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ABSTRACT 

One of the more significant error sources in using the 
for attitude determination of aircraft and earth orbiting satellites is multipath interference. 
The GPS signal is reflected fi'om the spacecraft structures resulting in a composite signal 
that differs fi’om the direct line-of-sight signal in terms of phase and amplitude. The 
resulting phase errors can have a significant impact on the attitude determination error. 
This thesis will describe the GPS signal generation and provide tools to. simulate the 
signal. A Kalman filter was developed to track the phase changes in the simulated GPS 
signals and its performance described in the absence and presence of multipath. 
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I. 


INTRODUCTION 


A. BACKGROUND 

In order for a spacecraft to perform its mission, it must be able to orient itself in 
relation to either the sun, the earth or both within an inertial reference frame. The primary 
system used to accomphsh this is the spacecraft attitude determination and control 
subsystem (ADCS). The ADCS uses external references to determine the vehicle’s 
absolute attitude with respect to some fixed inertial reference frame. Traditionally, these 
external references have included the Sun, the Earth’s infrared (IR) horizon, the local 
magnetic field direction, and the stars. Recent technological developments have added a 
new tool to the external reference arsenal - the Global Positioning System (GPS). 

Typical traditional ADCS sensor characteristics are shown below in Table 1. 
(Larson and Wertz) 


SENSOR 

TYPICAL PERFORMANCE 

RANGE 

WEIGHT (KG) 

POWER (W) 

Inertial measurement unit 

(gyros & accelerometers) 

Gyro drift rate = 0.003°/hr to 

l°/hr 

3 to 25 

10 to 200 

Sun Sensors 

Accuracy=0.005° to 3° 

0.5 to 2 

0 to 3 

Star Sensors (scanners & 

mappers) 

0.0003° to 0.01° 

3 to 7 

5 to 20 

Horizon Sensors 

0 

2 

0 

O 

2 to 5 

5 to 10 

Magnetometer 

0.5° to 3° 

0.6 to 1.2 

< 1 


Table 1 Typical ADCS Sensors 


In comparison, the Trimble TANS Vector attitude determination system using 
GPS is capable of an accuracy of 0.1° with a mass of 2.2 kg (GPS receiver and antennas) 
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and a power requirement of 7.5 W. Additionally, the Trimble system also provides three- 
dimensional orbital position information; which other sensors are incapable of providing. 

As described by Cohen (1996), attitude determination using GPS is achieved by 
determining the relative positions of multiple antennas, mounted to the vehicle, based on 
differences in the carrier phase measurements at the different antennas. The basic premise 
is that the GPS satellite is so distant, relative to the receiving antenna separations, that 
arriving wavefronts can be considered as effectively planar as shown below in Figure 1. 



Figure 1 Attitude Geometry 


The phase of a signal, ^, is a function of the wavelength, A,, and distance traveled. 


d: 



The wavelength is a fimction of the speed of light, c, and the frequency, f The 
distance traveled is a function of the speed of light and the time difference of arrival 
between antennas. At: 
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d = cA? 


Therefore a relationship between, time and phase can be shown to be 
cAt 



The incoming signal will arrive at the antenna closest to the transmitter before 
reaching the most distant antenna. By measuring the difference in carrier phase between 
the antennas, a receiver can determine the relative range between any two antennas. The 
measured differential phase, A<p, for baseline / and sateUite /, is proportional to the dot 
product of the baseline vector, x , and the line of sight unit vector to the transmitter, i. 
As shown in Figure 2, the GPS receiver measures only the fractional part of the 
differential phase. The integer component, k, must be resolved by independent means. 

The differential phase can therefore be expressed as: 



where ^ij is the noise in the measurement. Given the differential phase 


measurements between the different antennas, the spacecraft attitude can then be resolved 
as shown by Cohen (1995). Filtering this measurement will be the focus of this thesis. 

The dominant noise source in differential phase measurements is multipath, which 
occurs when the signal arriving at the antenna consists of the reflected signals in addition 


to the line-of-sight source. The reflected signal is phase shifted with respect to the original 


transmission and appears as additive noise at the antenna. Because the antenna locations 
are different, the multipath signature at each antenna is unique and will not be common to 
every antenna. Multipath accounts for more than 90% of the total error budget in carrier 
phase measurements. (Lightsey) 
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B. PURPOSE 


This thesis will discuss the generation of the GPS carrier signal which is used in 
making differential phase measurements. Models will be developed for simulating the 
multipath environment and the GPS signal. A Kalman filter which tracks the phase in the 
signal for use in the differential phase measurements wiU be described. Finally, results of 
the Kalman filter’s use in a multipath environment will be presented. 

Chapter n will discuss the GPS Signal Generation and present a model of the 
signal. The multipath interference will be examined and modeled in Chapter III. Chapter 
rv provides the derivation of the Kalman filter for tracking the phase of a GPS signal. 
Chapter V wiU show the results when the filter is used in tracking the GPS phase. The 
final chapter will summarize the findings and present other ideas for minimizing multipath 
and suggestions for future investigation. 
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n. GPS SIGNAL GENERATION 


The GPS signal is generated as shown in Figure 3. The GPS satellite transmits 
two carrier frequencies, LI (primary frequency) and L2 (secondary frequency). These 
carrier frequencies are modulated by spread spectrum codes with a pseudorandom noise 
(PRN) sequence unique to each satelUte. The LI frequency is modulated by the 
navigation message along with two PRN codes: the coarse/acquisition (C/A) code and the 
precision (P) code. The GPS L2 signal is modulated by one of the codes using Binary 
Phase Shift Keying (BPSK). The GPS LI signal is modulated using both codes by an 
unbalanced quadrature phase shift key (QPSK) modulator which is actually two BPSK 
modulators with different amplitudes. This development closely follows that given by 
Kaplan (1996). 

A. PN CODE GENERATION 

The PRN codes used in the generation of the GPS signal are non-maximum length 
codes also knovm as Gold codes. These codes are generated from an n-bit shift register 
generator. This register has an initial state or fill. After each clock cycle the output of a 
specified cell is used as the input to a modulo-2 adder (exclusive-or). The output of the 
modulo-2 adder is then fed back into a specified stage of the register. As an example, 
consider the shift register shown in Figure 4. This register uses the 3^** and 5* stages as 
the input to the modulo-2 adder and feeds the result back into the first stage. The 
polynomial used to describe the design of linear code generators is of the form 

i+Z^' 

where x' means the output of the i* cell of the shift register is used as the input of 
the modulo-2 adder and the 1 means the output of the adder is feed in to stage 1. The 
polynomial for the register in Figure 4 is then 1 + + X^. If this register was initially 

filled with all I’s at time 0, the resulting register for the first ten clock cycles would be as 
shown in Table 2. The binary signal output is taken from the last register and in this case 
would be 11111000110. 
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Information 


Figure 3 GPS LI andL2 Signal Generation 


Clock 



Output 


Shift Register 

Figure 4 PRN Code Shift Register for Example Polynomial 7+x'+x^ 
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Register 

time 1 2 3 4 5 

0 11111 

1 0 1111 

2 0 0 1 1 1 

3 0 0 0 1 1 

4 1 0 0 0 1 

5 110 0 0 

6 0 110 0 

7 10 110 

8 110 11 

9 1110 1 

10 0 1110 

Table 2 PRN Code Generation Example Results for Polynomial 7+x^+x^ 


The Matlab m-function gen poly.m can be used to generate the PRN sequence for 
any user defined polynomial and register. 

B. C/A CODE GENERATION 

The GPS coarse/acquisition (C/A)-code is a Gold code generated by two 10-bit 
shift registers, Gl and G2. The polynomial for each is given below (Kaplan, 1996) 

G1 = 1 + + X’® 

G2 = 1 + X^ + X^ + X® + X* + X® + X*® 

Both C/A-code registers have an initial fill of all 1 ’s. Each code is derived and the 
result of G2 is then delayed by the satelUte PRN code number. This delay effect is 
obtained by the exclusive-or of the selected positions of the G2 register. The C/A-code is 
then obtained by the exclusive-or of the Gl and delayed G2 signals. An example for 
satellite vehicle 1 is shown in Figure 5. 

The code tap positions and initial code sequences are given in Table 3. The m-files 
gl.m and g2.m generate the codes for a specified satellite vehicle number. 
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C/ACode 

Gl(t) 


Figure 5 C/A-Code Generator for Satellite Vehicle 1 with C/A Code Tap Selection 206 

(from Kaplan) 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 


C/A-Code 

Tap 

Selection 

C/A-Code 

Delay 

(Chips) 

P-Code 

Delay 

(Chips) 

First 10 C/A 

Chips 

(Octal) 

First 12 P- 

Chips 

(Octal) 

2©6 

5 

1 

1440 

4444 

3©7 

6 

2 

1620 

4000 

408 

7 

3 

1710 

4222 

509 

8 

4 

1744 

4333 

109 

17 

5 

1133 

4377 

2010 

18 

6 

1455 

4355 

108 

139 

7 

1131 

4344 

209 

140 

8 

1454 

4340 

3010 

141 

9 

1626 

4342 

203 

251 

10 

1504 

4343 

304 

252 

11 

1642 

4343 

506 

254 

12 

1750 

4343 

607 

255 

13 

1764 

4343 

708 

256 

14 

1772 

4343 

809 

257 

15 

1775 

4343 

9010 

258 

16 

1776 

4343 

104 

469 

17 

1156 

4343 

205 

470 

18 

1467 

4343 

306 

471 

19 

1633 

4343 

407 

472 

20 

1715 

4343 

508 

473 

21 

1746 

4343 

609 

474 

22 

1764 

4343 

103 

509 

23 

1063 

4343 

406 

512 

24 

1706 

4343 

507 

513 

25 

1743 

4343 

608 

514 

26 

1761 

4343 

709 

515 

27 

1770 

4343 

8010 

516 

28 

1774 

4343 

106 

859 

29 

1127 

4343 

307 

860 

30 

1453 

4343 

308 

861 

31 

1625 

4343 

409 

862 

32 

1712 

4343 


Table 3 GPS Code Phase Assignments (from GPS-ICD-2000) 
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test CA.m is a script file to test C/A code generation. This program calculates 
the first 10 C/A-chips and converts the result to octal numbers that match the First 10 
C/A-Chips (Octal) of the GPS-ICD-2000 table shown in Table 3 

For example the first 10 chips of the PRN 5 C/A-code is 1001011011 which is 
1133 in octal. 

C. P-CODE GENERATION 


The GPS Precision (P) Code is generated by PRN sequences using four 12 bit shift 
registers whose polynomials and initial states are shown below: 


Register 

Polynomial 

Initial State 

XIA 

1 +x*+x^+x^+x* + x^+x^®+x" + x*^ 

001001001000 

XIB 

1 +x* + x^+x^ + x^+x^ + x’+x*+x^+x^''+x“ + x‘^ 

010101010100 

X2A 

1 +x*+x^+X‘+x^+x’+x*+x^+x*®+x"+x*^ 

100100100101 

X2B 

1 + X^ + X^ + X'* + X* + x^ + x*^ 

010101010100 


Table 4 P-Code Generator Polynomials and Initial States 


The P-code is 7 days in length at a chipping rate of 10.23 MBps. The 7 day 
sequence is the Modulo-2 sum of two subsequences XI and X2i; their lengths are 
15,345,000 chips and 15,345,037 chips respectively. 

XI itself is generated by the Modulo-2 sum of the output of two 12-stage registers 
(XIA and XIB) short cycled to 4092 and 4093 chips respectively. When the XIA short 
cycles are counted to 3750 the XI epoch is generated. This occurs every 1.5 seconds, 
after 15,345,000 chips of the XI pattern have been generated. The XI period is defined 
as 3750 XIA cycles (15,345,000 chips) which is not an integer number of XIB cycles. To 
accommodate this the XIB shift register is held in the final state (chip 4093) of its 3749*** 
cycle for 343 chips. The X2 sequence is generated in a similar manner, however the X2 
period is defined to be 37 chips longer than the XI period in order to force the X2A and 
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X2B epochs to precess with respect to the XIa and XIB epochs. This precession 
continues until the end of the 7 day period. 

The Y-code’ is used when the anti-spoofing mode of operation is activated. It is 
generated in the same fashion as the P-code, however the polynomial generating equations 
are classified. (ICD-GPS-2000) 

xla.m, xlh.m, x2a.m and x2h.m are used to calculate the PN codes. The m-file 
test P.m is used to test P-code generation. 

D. DATA GENERATOR 

The output of the data generator is the navigation message, D(t). The message 
includes satellite vehicle ephemerides, system time, clock bias information, status 
messages, and C/A to P-code handover information. The navigation message is generated 
at 50 bps and is added to the P and C/A-codes to modulate both the LI and L2 signals. 
aCD-GPS-200) 

E. BPSK MODULATOR 

After the various codes are generated they are then used to modulate the carrier 
signal using binary phase shift keying, BPSK, resulting in the LI or L2 signal. The 
Simulink program used to simulate the BPSK modulator is shown below in Figure 6. 


' The Y-code and P-code are sometimes discussed together and are referred to as P(Y)-code meaning P (or 
Y) code. 
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Inport 


►JV—® 


Zero Order Hold 


Relay 




Sine Wave 


Product 


Mux 


> t BPSK 


To Workspace 


Mux 


Outport 


Figure 6 BPSK Simulink Diagram 

The PRN code is input via the inport. It is sent through a zero-order-hold and 
relay function to generate a non-return to zero (NRZ) signal, which is multiplied by the 
sine wave, resulting in a binary phase shifted signal. 
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BPSK Simulation 
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Figure 7 BPSK Simulation Results 


Figure 7 shows an example of the BPSK simulation results. This figure shows the 
results: NRZ signal in the top graph, the sine wave in the second graph, the product of the 
two, or the BPSK signal in the third graph, and a composite of the data and resulting 
signal in the bottom graph. 


F. LI SIGNAL 

The P-code at, 10.23 x 10® chips per second, and C/A-code, at 1.023 x 10® chips 
per second, are combined with the 50 bits per second data and are then used to modulate 
the LI frequency (154 • 4 = 1575.42 MHz ) as shown in Figure 3. This results in an 
Unbalanced Quadri-Phase Shift Keying (UQPSK) modulation with the data bits added to 
the ranging codes, the C/A-code signal lagging the P-code signal by 90°; and the C/A- 
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code signal power nominally exceeding the P-code signal power by 3 dB. (Struza) 
Mathematically, this can be represented by: 

= ^j[P(/)® D(/)]cos((y,0 + V2/l,[G(/)®/)(0]sin(®,t) 
where P(t) is the P-code, D{t) is the data, and G(t) is the C/A-code. 

G. L2 SIGNAL 

The L2 frequency (120- = \221.6MHz) is modulated by either the exclusive-or 

of the P-code and data, the P-code, or the exclusive-or of the C/A code and data, as 
selected by the control segment. Mathematically, this can be represented by: 

4 (^20 = A [^(0 ® Z)(0] cos(ty20 
4(^20 = 4[^(0]cos(fi)2/) 

4(<y20=4[^(0 ]cos(<»2/) 


H. SUMMARY 

The GPS signal consists of two separate signals, LI and L2. These signals are 
generated using BPSK and unbalanced QPSK modulation. The modulation signal is 
generated using the P and C/A codes to which the data signal is added. The Matlab files 
init.m, sig gen.m and my glb.m are used to initialize and generate a sample of the codes 
used to modulate the signal. These codes are then used by the Simuhnk programs 
LI sim m and L2 sim .m to generate a sample of the signal. In the next chapter the 
multipath environment will be described. 
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m. MULTIPATH 


A. INTRODUCTION 

This chapter will discuss the multipath environment that impacts the use of GPS 
for attitude determination and will present a method to model multipath interference using 
Matlab. 

Multipath occurs when the signal arrives at the antenna from reflected surfaces and 
from the line of sight sources. A simple case with a single reflected path signal is shown in 
Figure 8. 



The reflected signal is phase shifted with respect to the original transmission and 
appears as additive noise at the antenna. When antenna locations are different, the 
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multipath signature at each antenna is unique and the error in the phase measurement is 
not common. (Lightsey) Multipath interference, due to environmental factors, will impact 
the signal, however this effect will be common to each received signal. 

For the developed model it will be assumed that the multipath environment will be 
static. In reality the multipath environment will be varying due to the satellites attitude 
and position relative to the GPS satellites changing over time. 

B. MATHEMATICAL MODEL 

In the absence of the multipath, the input signal to a GPS receiver for the LI 
QPSK signal is given by: 

s(t) = cos(a>t +a(t)f) 

where is the received signal amplitude and a(t) = 0,1,2 or 3 is the pseudo¬ 
random code waveform that phase modulates the carrier. (Kumar) In the presence of N 
multipaths, in addition to the direct line of sight path, the input signal is given by 

^m(0 = ccAi) + - r,)+.. .+a^s(t - ) 

where a. and t, denote the amplitude and delay of the i* multipath for 
i=l,2,.. .,N. Graphically, this can be represented as shown in Figure 9. 
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Figure 9 Multipath Block Diagram 
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Typically for BPSK and QPSK signals the carrier frequency operates at a much 
higher frequency than the modulation rate. As a result, a small value of delay, t , vrill 
result in large phase differences. When the multipath signal with a large phase difference 
is added to the direct path there may be destructive adding which can significantly impact 
the received signal. 

C. SIMULATION MODEL 

In the model mpathm, three vectors are used to develop the signal 

delay -1 by vector of delayed samples. The first entry is the undelayed 
signal, the next entry is the signal delayed by one sample period. The length of this vector 
is the maximum value in the tau vector and is initialized as all zeros. 

tau - vector of delay times. Each entry is 1 minus the delay time in sample 
periods, dt, and is used to index the delay vector. The first entry is 1 and represents a 
delay of zero. 

alpha - vector of amplitudes. Each entry is the amplitude of the 
corresponding delayed signal. The first entry is a 1 in order to return the undelayed hne 
of sight signal. 

As an example, the BPSK signal shown in Figure 10, with N=3 multipath signals, 
will be used to demonstrate the use of the vectors and the algorithm used to develop the 
signal 

- ■Z’o) + - Ti) + ^2 - ^ 2 ) + 

where 

T = [l 10 12 20 ] 
a = [\ 0.5 0.3 0 . 1 ] 

thus 

(0 = -^(0 + 0-5 's(t-9' dt) + 0.3 • -11 • + 0.1 • -19 • dt) 
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If the sample period, dt, equals 0.01 the following signal is produced using the m- 
file mpath.m 


Multipath simulation 
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Figure 10 Multipath Example Signal 

The top graph in Figure 10 shows the undelayed signal, the middle graph shows 
the direct path and multipath signals plotted individually and the bottom graph shows the 
sum of the individual and multipath signals. 

If the multipath delay time and amplitude are increased with the following 
parameters 
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combined composite undelayed 


T = [l 20 22 30] 
a = [l 0.5 0.3 0.5] 

the resulting multipath signal will be as shown here; 
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Figure 11 Multipath Example with Large Delay 

Notice that in both examples the multipath signal has a different amplitude from 
the original signal and that the larger the delay the larger the change in the phase of the 
signal. Below is an example of destructive interference obtained with the values 
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combined composite undelayed 


r = [l 50 52] 
a = [l 0.5 0.4] 
Multipath simulation 



Figure 12 Multipath Example with Destructive Interference 

This chapter has described the multipath environment and presented a simulation 
model to generate a signal based on user defined parameters. The next chapter will 
describe a Kalman filter which is used for tracking the GPS signal. 
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rV. KALMAN FILTER FOR TRACKING PHASE 


A. DISCRETE KALMAN ESTIMATOR 

The development of a filter to track the GPS signal in a multipath environment 
began with the implementation of a Discrete Kalman Estimator for tracking the phase 
changes of BPSK signal. The derivation below closely follows that given by Professor 
Kirk (1975). 


1. Plant Equations 


The plant equations were derived as follows: 

X, = COS(fiJ t) 


Xj = X, = -a sin(<u t) 


x(0 = 






^ cos(o)t) ^ 
V-n; sin(6; r)> 


i(0 = 




\X2J 


^ -n;sin(<ur) 


^ 0 r 

^ COS(CL)t) ^ 

V-<y^ cos(n;0> 


K-co^ 0; 

V-ctJ sin(fy t)J 


= Ax(0 


2. Plant Characterization 


The discrete model of the plant is characterized by the linear discrete difference 
equations: 


x(k +1) = ^ x(^) + w(k) 


z(k) = Cx(k) + v(k) 

where: 

x(k) is the n-dimensional state vector at time k, 

z(k) is a q-dimensional output measurement vector at time k 

w(k) is a p-dimensional vector of random forcing inputs at time k 
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<}> is the state transition matrix 
C is the observation matrix 
and 

v(k) is the q-dimensional vector of random measurement noise at time k 

3. Gain, and Covariance Equations 

The gain, G, and covariance equations derived by Kirk (1975) are: 

G{k) = P{k\k -1) • C"” • [C • P{k\ k-\)C^ + 7?]'' 

P(Jc\k) = [/ - G(yt) • C] • P{k\k - 1) 

P{k + \\k) = <l)-P{k\k)-<1)^+0 

where 

C is the observation matrix 
I is the identity matrix, 

G is the Kalman gain 
P is the variance of estimation error 
Q is random process forcing input covariance matrix 
R is the variance of measurement noise. 

The estimator equations are: x{k\k) = x(k\k -1) + G(k)■ [z(k) - C• x(k\k - 1)] ^ 
x(Jc\k-\) = (j>-x(Jk - 1|A: -1) + A- w(A:-1) 

with the initial condition x(0|-l) = x^ 


^ ^ given j and is defined as the estimate of x at time j given measurements up to and 

including time j. 
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4. Assumptions 


a) The measurement noise has zero mean; 

£[v(i)]=0 /:=0,1,..= 

is uncorrelated and has covariance R(k): 

£[v(i)v'0)] = £[vOV(i)] = 
where 5^^ is the Kronecker delta function defined by 

[0 j^k 

b) The initial state is a random variable with known mean 

£[)t(o)]=;; 

and covariance; 

^[[^(0) - [x(0) - ] = M 

The measurement noise and initial state are uncorrelated; 

£[x(0)v^ (iS:)] = E\v{k)x^ (0)] =0 k = 0,1,... 

c) by the linear relationship 

x{k\k) = x{k\k -1)+ Cj{k) ■ [z{k) -C-x{k\k- 1)] 

where G(k) is the Kalman gain and z{k)- C • x{k\k — l) is the correction term, 

the difference between the measured and predicted value. 

d) The random process forcing input has zero mean 

£[w(A:)] = 0 k = 0,\,... 

is uncorrelated and has covariance Q(k) 

£[w(*)«."0)] = £[wO)M-’'(i)] = Q{k)S^ J,k = 0,1,2..,. 

e) The forcing random process and initial state value are uncorrelated 

^E[] is (he expected value or mean of the random variable. 
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£'j^w(A:)x^(0)j = £j^Jf(0)>v^(A^)j = 0 = 0,1,2,... 

f) The forcing random process and measurement noise process are uncorrelated 
£[w(4')v"0)] = £[vO)w"(i)] = 0 t = O.U,... 

COMPUTATIONAL PROCEDURE 

The computational procedure is shown in the following algorithm. 

Initialize with 
k=0 

x(0|-l) = ^ 

P{k\k-\) = M 

The procedure iterates on k in the loop; 

G{k) = P(A:| A: - 1 ) • • [c • P{k\ A: - 1 ) • + rJ' 

P(k\k) = [/- G(k) ■ C]P{k\k - 1) 

P{k + 1|A:) = P(A:|A:) -^^+0 (this will be P(A:|Ar -1) next time through) 

take measurement z(t) 

x(k\k) = x(A:| A: -1) + G(k) ■ [z(k) - C ■ x(k\k- 1)] 
jc(A^|A:-l) = ^•x(A^-1|A^-1) + A-m(A:- 1) 
k=k+l 

End of the Loop 


PHASE TERMS 


The QPSK signal can be represented by 


s(t) = Acos((i)t + <^{t)) 

where ^(/) represents the phase change of the signal 


^(0 = 


<Po 

<l>o+- 
° 2 

4>^ + n 

, "in 

+ - 

° 2 


Substituting in these phase values and simphfying using trigonometric identities 
results in 



Acos{(ot + ^g) 
-Ai\n{(ot + 
-Acos,{(ot + 
A%m{(ot + ^g) 


Using these relationships and the sine and cosine information in the state estimates, 
the Kalman filter can estimate which phase change has occurred. The error term is 
calculated by fiaiding the estimate of the signal that minimizes the difference between the 
measured value. Assuming 

^ C0S{0}t + <f>g) 


x(A:|A--l) = 


-(osm{cot + (j>„)) 


C = [1,0] 
C,., = [0,1] 


the four error terms can be calculated as shown below, and the minimum value is 
used to update the prediction 

= z{k) - C'L(Jk\k -1) = z{k) - cos(<y k + (j)^) 
e, = z{k) + C^^x(^| A: -\)l a = z(k) - sin(ru k + <l>^) 

= z(A) + Cx(A|A - 1) = z(A) + cos(m k + <pg) 

= z(A) - C^^x(A| A -\)l a> = z{k) + sin(<u A + 
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These error terms are then used to update the estimate of the signal without any 
phase shifting which will be called y, and an estimate of the phase shifted signal called 
y_ps. 

Next, the results of using the developed Kalman filter to track the phase of a 
BPSK, QPSK, and unbalanced QPSK signal will be shown. 
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V. 


INITIAL RESULTS 


This chapter will show the results obtained when the Kalman filter, discussed in the 
previous chapter, is used the track the phase of a BPSK, QPSK and unbalanced QPSK"* 
signal. The performance of the Kalman filter, will be shown with and without noise, and 
in the presence of multipath interference. It will be shown that the Kalman filter will track 
the combined line of sight and direct path signals phase and amplitude. It will then be 
shown that when given a priori information the filter can distinguish between the line of 
sight and multipath signals and correctly track them, given a static multipath environment. 

A, SIGNAL GENERATION 

The Matlab m-file uijyars.m sets up the user interface shown in Figure 13 that can 
be used to control the simulation. 

Each entry can be edited by the user to change each of the parameters discussed 

below; 

• Signal freq - Sets the fi’equency, of the oscillator, in Hertz, used by the 
carriers in the BPSK and QPSK signal generation. 

• Phase Offset - Determines the initial phase offset, of the line of sight signal in 
degrees. 

• Amplitude - Sets the amplitude of the generated signal. 

• - This sets the sampling period in seconds. In the example below Ts =0.01, 
resulting in 100 samples per second. 

• Stop time - Sets the Simulink parameter stop time to control length of 
simulation. 

• Tb- Sets the bit period which determines the length of time for zero-order hold 
per data bit. 


'' An unbalanced QPSK signal is a QPSK signal with differing energy in the in-phase and quadrature 

phase components. For GPS the LI signal the amplitude of the quadrature term is equal to V2 times the 
in phase term. 
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SNR dB - Using this value and the entered amplitude value of the signal, the 
random number generator is multiplied by a scaling factor to simulate the 
desired SNR level. 

tau - Delay vector used to control the multipath delay times. See Chapter III 
Section C for more information. 

alpha - Amplitude vector - Used to control the multipath amplitudes. 



Figure 13 User Interface to Control BPSK and QPSK Simulation 






The pulldown menu Init initializes all the entered values. Plot gives the user a 
choice of plots to choose from to document the output of the simulation. Sim allows the 
user to choose which simulation will be run. 

B. BPSK 

Using Matlab and Simulink, a BPSK signal was generated and the algorithm 
discussed in the previous chapter was used to track the phase changes in the signal.^ 

1. BPSK Signal Without Noise 

Figure 14 shows the performance of the algorithm when there is no noise present. 
In this graph, x is the plant state, z is the measured signal (x +noise), y is the recovered 
carrier, and y_ps is the recovered phase shifted signal. The last graph shows the error in 
the recovered signal which is the difference in the generated signal and recovered estimate, 
x-y_ps. Notice there is an initial error during initialization, however the error quickly goes 
to zero as expected in the absence of noise. 

2. With Noise 

Keeping all variables the same as above and adding zero mean Gaussian white 
noise, by using the random function in Matlab, with an SNR value of 12 dB®, the carrier is 
still recovered, however small errors are now introduced in the estimate of the phase 
shifted signal. 


^ See Matlab m-files kaljnit.m, kalmanl.m and the Simulink program kaljsim 
® Spilker states that typical values for C/No are 40.6 dB-Hz, wldch corresponds to an SNR of 25 dB. See 
Appendix C for a link budget estimate of the signal to noise ratio. 
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Figure 14 Kalman Filter as Phase Lock Loop for BPSK Without Noise 


34 

















Kalman Filter 



time ~ sec 

Figure 15 Kalman Filter for BPSK with Noise 


C. QPSK 

Next the results for a QPSK simulation are shown. These plots were generated 
using kalqsim.m 

1. Without Noise 

Figure 16 shows the performance of the algorithm when there is no noise present. 
The true signal, x, and the estimate y_ps and the error in the estimate, x-y_ps are all 
plotted. 
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Figure 16 Kalman Filter for QPSK, Without Noise 


2. With Noise 

Figure 17 shows the estimation when white noise^ is added with an SNR value of 
12 dB. Figure 18 shows the error in the estimation. The carrier is still recovered, however 
small errors are now introduced in the estimate of the phase shifted signal. The cause of 
the error will be explained in the next section. 


’ Throughout this study only gaussian white noise was considered. Colored noise from other sources, such 
as electromagnetic interference was not considered. It was assumed that any such noise would have the 
same affect at both receivers and would not change the difference in the measured phase values. 
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Kalman Filter 



Figure 1 7 Kalman Filter for OPSK with Noise 


Kalman Filter Error 




Figure 18 Error in OPSK Estimator 
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D. SOURCES OF ERROR IN NOISY CASE 


The error in the estimation caused by the noise can best be explained by examining 
the BPSK simulation. In the algorithm, the Kalman filter tracks the carrier signal and 
estimates which generated phase shifted version has the minimum error when compared 
with the measured signal z. Figure 19 is a plot of the original BPSK signal, x, the signal 
with noise added, z; the estimated signal, y_ps; the two possible signals, xO and x2, that 
must be decided between; and the error, eO and e2, used as the basis for the decision. 

The error in the signal is best described as “clicks” which occur when the Kalman 
filter is unable to distinguish between the correct and incorrect phase because they both 
have approximately the same value. The reader is invited to examine this effect between 
times 1 and 2 seconds. The correct estimated signal is X2 which results in the minimum 
error over the bit period, however at approximately time 1.9 the out of phase signal 
estimate XO has the same value as X2 and the error eO is less than el due to the added 
noise. This produces an erroneous decision to switch phase. Typically the erroneous 
decisions all occur when both the correct and incorrect signal estimates are near a zero 
crossing. This is due to the nature of the cosine fimction which results in 
cos(90°)=cos(270°)=0. 
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Figure 19 Kalman Filter Error Terms and Decision Options 
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E. UNBALANCED QPSK 


As discussed in Chapter II, the GPS LI signal is actually an unbalanced QPSK 
signal with differing energy in the quadrature and in-phase components. As shown in 
Chapter II Section F on page 15 the LI signal is given by 

Lj (<y,0 = [P(0 ® D(0] cos((yit) + V2y4, [G(/) <8> D(/)] sin(cy,/) 

An example of the unbalanced QPSK signal is shown in Figure 20. The top graph 
shows the composite of the binary data representing [G(/) 0 D(0] with the quadrature 

component, the middle plot shows the binary data representing [P(O0Z)(Q] with the in- 

phase component. In both cases the Ai has been normalized to 1. The last plot shows the 
sum of the two signals which is the resulting unbalanced QPSK signal. 


QPSK Simulation 



Figure 20 Unbalanced QPSK Signal Generation 
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When the Kalman filter is used to track this signal as shown in Figure 21, errors 
occur when the in-phase and quadrature terms have differing data. The filter locks on to 
the initial amphtude and assumes the amplitude will not change. In this example the 
absolute value of the amplitude is the same when the data bits are the same (both +1 or 
both -1) and is different with the data bit have opposite values. The state estimation was 
developed under the assumption the signal would have a constant amplitude. In order to 
resolve this problem the logic for phase change detection would need to be modified. 
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Figure 21 UQPSK Without Noise 
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When noise is added to the unbalanced QPSK as shown in Figure 22 the 
performance is the same as in the QPSK example when the data bit for the in-phase and 
quadrature phase term is the same, however when they are different, the performance 
degrades. The error in the estimation is shown in Figure 23. 
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Figure 22 UQPSK with Noise 
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Kalman Filter Error 




Figure 23 Error in UQPSK Estimator 


F. MULTIPATH RESULTS 

1. Without iVio/7 Knowledge 

The multipath parameters used to generate Figure 10 are used to generate the 
measured signal z and the results obtained are shown below in Figure 24. The estimate, 
y_ps, and the direct path signal, x, are shown plotted together and the error, or the 
difference between x and y_ps is also shown. 
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time ~ sec 


Figure 24 Results of Kalman Filter Estimation in Presence ofMultipath 


The Kalman filter has correctly determined the phase changes, however due to the 
multipath, the amplitude of the estimation is now different fi'om the original signal. Using 
the larger multipath spread parameters used to generate Figure 11 the resulting phase 
error is observed as shown Figure 25. 
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Kalman Filter Error 



Figure 25 Results of Kalman Filter Estimation in Presence of Multipath with Increased 

Delay 


The results of these last two plots indicate the Kalman filter performance in the 
presence of multipath is much worse than in the presence of white noise. If multipath is 
present it will introduce a phase error that will be unique to each antenna. When the 
differential phase is calculated it will now include the multipath phase error and will no 
longer be just a function of the geometry. 
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2 . 


With A Priori Knowledge 


If the multipath amplitude and time delay are known, the Kalman filter can be 
modified to track both the the direct path signal and delayed multipath signals. A simple 
case with a single multipath channel is shown below with the following assumptions; 

The combined signal is given by: 

The incoming signal is not changing phase and is given by 

5(0 = A cos(o) t + ${t)) 

where 0(t)= a{t) f and a(t) = 0,1,2,or 3 representing the intelligence for the BPSK 
or QPSK modulation. 

The time delay r is converted into the phase delay <j> 

s{t-x) = Acos{(o{t-r) + d{t-t)) = A cos{a)t - o)t += Aco%{o)t + ^) 
and the composite multipath signal is now given by 

5 „ (0 = ^ cos(<u t + d) + a A cos(ty / + ^ ) 

The plant equations are now 

X, = A cos(<y t + &) 

Xj = Xj = - 0 ) A sin(fy t + 6) 
x^= a A cos(<y t + (j>) 

X4 = Xj = -a a A sin(o t + (j>) 


x(0 = 




^ A cos(o t + 6) ^ 

Xj 


-CO A %\xi{cot + d) 

^3 

— 

a A cos(<y t + ^) 

U 4 ; 


G)a As\n{(ot + (l>)) 




^ -G> A sin((D t + d) ^ 


^ 0 1 0 0^ 


^ A cos(iO / + ff) '' 



-co^Aco%{o}tO) 


-co^ 0 0 0 


- caA sin(a> t + 0 

^3 


-CO a A sin(6}/ + ^ ) 


0 0 0 1 


a A cos{co t + ^) 

U4J 


K-co^a A cos(<y t 


\ 0 0 -6)^ 0> 


CO a A siniot + ^)) 


= Ax(0 
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Since a priori knowledge is known the values for the estimate at time zero can 
now be initialized as follows: 


x(Ol-l) = 


A 

0 

a A cos{<f>) 
-coa A sm{(p) 


Figure 26 shows the line of sight signal, delayed signal, and the measured signal z 
which were generated using a=0.5 and t=10 with Ts=. 1. For this example the multipath 
environment is assumed to be static. In reality, the multipath environment would be much 
more dynamic due to the satellites orbit, the satellites changing attitude, and the changing 
position of the GPS satellites. 


to 



N 



Figure 26 Example Signal for A Priori Case 
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los & xkkps(1,;) (estimate = 



Figure 27 A Priori Knowledge Results Showing Signals and Estimates 


Figure 27 shows the estimations and original signals plotted on top of each other. 
This demonstrates how, given a priori knowledge, the line of sight and delayed signal can 
be tracked. Now there is no phase error due to the multipath signal, unlike the examples 
shown without a priori knowledge when the filter tracked the phase of the combined 
signal. 

When noise is added in this simple case with 12 dB SNR as shown in Figure 28 
and Figure 29 there is also improved performance with a priori knowledge. 
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Figure 28 Generated Signal for A Priori Example with Noise 


los & xkkps(1.;) (estimate = 



Figure 29 A Priori Knowledge with noise Residts 
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G. EFFECT OF Q MATRIX ON KALMAN GAINS 

The random process forcing input covariance matrix Q is a measure of the 
expected excitation of the noise. It is used in the calculation of the Kalman gains as 
shown in Chapter IV, Section A-3 on page 26 which are restated here: 

G{k) = P{k\k-\) C^ [C-P(J(\k-\)C^+ i?]'' 

P(k\k) = [/- G(k) • C\P{k\k - 1) 

P{k-¥\\k) = (j>-P{k\k)-<j>^+Q 


The Q matrix is a function of the state noise coefiBcient matrix, F, and the 
variance of the forcing function, w, as shown below. 


df 



dt 


Q = T •var[>v] 


dt'^ 



2 


dt^ 

dt^ 


• var[w] 


As there is no forcing function, w, the value of Q was initially set to zero, however 
in order to determine its impact, different values for the variance of w, var[w], were used. 
The filter results and Kalman gain values G1 and G2 are shown in Figure 30 through 
Figure 33 for different values of var[w] and Q. As the value for var[w] is increased, the 
steady state gain increases. The only noticeable effect occurs when var[w] is greater than 
100. With an increase in var[w] and Q, there is also an increase in the Kalman gains. This 
results in more responsiveness. If the gain values are too large the filter will be too 
responsive resulting in more errors. 
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Kalman Filter - Signal Options 



Figure 30 


Ejfect of Gain Values with var[w]=0, 


Q = 


0 0 
0 0 


51 





















Kalman Filter - Signal Options 



Figure 31 


Effect of Gain Values with var[w]=lO, 0 = 


0.0003 

0.005 


0.005 

0.1 
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Kalman Filter - Signal Options 



Figure 32 


Effect of Gain Values with var[w]=100, 0 = 


0.0025 0.05 
0.05 1 
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Kalman Filter - Signal Options 



Figure 33 


Effect of Gain Values with var[w]=150, Q = 


0.0038 0.075 
0.075 1.5 
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Steady State Gam Values 



Figure 34 Steady State Kalman Gains mth varfwJ=J50 


The steady state gain values after time t=2 seconds are shown in Figure 34. For all 
other values of var[w] shown, the corresponding steady state Gain values approached 
zero. 

H. SUMMARY 

This chapter has presented the results of using the Kalman filter to estimate the 
phase in BPSK, QPSK, and Unbalanced QPSK signals to simulate the various formats of a 
GPS signal. Results have been shown for the performance with and without Gaussian 
white noise and with multipath interference. 

The designed filter works best with the BPSK signal. Gaussian white noise 
degrades performance due to the effects of the correct and incorrect phase having the 
same value at a zero crossing. It may be possible to modify the phase switching logic to 
reduce this effect by allowing a phase change to only occur at or near the expected time of 
the end of a bit period or when the error differences are greater than some minimum value. 
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The multipath interference was shown to result in errors in the estimate’s phase 
and amplitude. This error is caused by the filter tracking the phase of the combined 
multipath signal. It was shown that with a priori knowledge of multipath delay and 
amplitude, the filter is able to discriminate between the direct path and the delayed signals 
and to eliminate this error. Sources of the a priori information might include ground 
testing by placing a receiver into a self-survey mode and allowing the GPS satellites to 
pass overhead. (Lightsey) It may also be possible to use the technique suggested by 
Kumar and Lau to estimate the multipath parameters through seismic deconvolution 

For this case it was assumed that the multipath environment was not changing. In 
reality the multipath environment would be very dynamic as a result of the satellites orbit 
and changing attitude. Additionally, the attitude determination and control system only 
uses 4 out of possibly 11 visible satellites. The satellites selected would also effect the 
multipath environment. To overcome this may require fi’equent re-initialization of the 
Kalman filter. 
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VI. CONCLUSION 


A. SUMMARY 

This study has focused on the carrier phase signal which is used by GPS based 
attitude determination systems. A GPS carrier signal and a multipath environment were 
described and modeled. A Kalman filter was developed to track the GPS signal and 
results were shown for Gaussian noise and in the presence of multipath. 

The developed Kalman filter was shown to work best for the L2, BPSK signal due 
to the minimum complexity involved in the phase changing logic. It was shown that the 
filter is capable of discriminating between the line of sight and delayed signals and 
performed well in the presence of multipath, given good a priori knowledge. 

More refinement is required in the phase change detection logic to eliminate errors 
caused by the error terms having the same value at a zero crossing. The Q excitation 
matrix needs considerable adjustment to handle changing phase in the presence of 
multipath. The source of the a priori information needs to be developed. 

B. OTHER METHODS OF MULTIPATH MITIGATION FOR POSSIBLE 

STUDY 

As multipath accounts for 90% of the error in an attitude determination system, 
much work is currently being done in this area. Three categories of multipath mitigation 
methods that have been used include: 1) Methods involving the geometry of the 
multipath, 2) Analyzing the consistence of the different code measurements and 3) 
Exploiting information coming from different channels. 

Of these three, the easiest to implement involves geometry. One approach is to 
focus on improving the antenna or its position. This may be accomplished by placing the 
antenna in isolated locations minimizing the multipath, something that is not always 
possible on a satellite. Another method of multipath rejection is achieved by reducing the 
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left-hand circularly polarized gain (LHCP) of the antenna without reducing the right-hand 
circularly polarized gain (RHCP). This is because the satellite signals are RHCP and the 
reflected multipath signals tend to be either linearly polarized or even LHCP, depending 
on the reflecting surface. (Van Dierendonck) An example of an antenna designed 
specifically to reduce multipath effects is the choke ring antenna. This antenna is 
comprised of vertical-aligned concentric rings centered about the antenna element that are 
connected to the ground plane. The vertical rings shape the antenna pattern such that the 
multipath signals incident on the antenna at the horizon and negative elevation angles are 
attenuated (Kaplan). 

Other methods of dealing with the multipath problem involve modification to the 
GPS receiver for code measurements. One approach involves reducing the early-late 
delay spacing among the correlators in the GPS receiver code lock loop. This technique 
provides limited results if the early-late spacing is smaller than the initial delay error due to 
multipath. (Kumar) 

An example of the third category is the method presented by Kumar and Lau 
involving the use of optimal deconvolution to estimate the impulse response of the 
effective multipath chaimel and obtain an inverse filter which equalizes the multipath 
channel (Kumar) Another technique involves adaptively estimating the spectral 
parameters (frequency, amplitude, phase offset) of multipath in the associated SNR and 
then constructing a profile of the multipath error in the carrier phase. A multipath 
correction is then made by subtracting the profile from the actual phase measurement. 
(Comp) 
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APPENDIX A 


MATLAB PROGRAMS 


A. DO KAL.M 


% File Name: do_kal.m 
% Last Updated: 18 Apr 97 
% Written B>’: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to set up and nm kalman filter simulations 

clear 

uivars 

% End of File do_kal.m 


B. G1.M 


% File Name: gl.m 
% Last Updated: 18 Oct 96 
% Written B)’: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the Gl(t) signal 
% for GPS C/A Code generator 
% 

% 


function [new_R, g] = gl(R, xl) 


if xi==l 

R = ones(size(R)); % Reset to all ones at XI epoch 
end; 

new_R=gen_poly([l,3,10],R); 
g=R(10); 


% End of File gl.m 
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c. 


G2.M 


% File Name: g2.m 
% Last Updated: 18 Oct 96 
% Written By; T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description; This function generates the G2(t) signal 
% for GPS C/A Code generator 
% 

% 


flmction [new R, g] = g2(R, xl, sv) 


ifxl==l 

R = ones(size(R)); % Reset to all ones at XI epoch 
end; 

new_R = gen_poly([l 2 3 6 8 9 10], R); 

%new_R(l) = mod2(R(2)+R(3)+R(6)+R(8)+R(9)+R(10)); 
%new_R(2:10) = R(l:9); 


% phase select logic, see p. 92 KAPLAN 
% This changes based on satellite choosen. 
sv_tap =[26 
37 
48 
59 

1 9 

2 10 
1 8 
29 
3 10 
23 
34 
56 
67 
78 

8 9; 9 10;1 4;2 5;3 6; 4 7; 5 8; 6 9; 1 3; 4 6; 5 7; 6 8; 7 9; 
8 10; 1 6; 2 7; 3 8; 4 9; 5 10; 4 10; 1 7; 2 8; 4 10]; 


g=xor( R(sv_tap(sv,l)), R(sv_tap(sv,2))); 


% End of File g2.m 
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D. 


GEN POLY.M 


% File Name; gen_poly.m 
% Last Updated; 18 Oct 96 
% Written By; T. H. Newman 
% Operating Environment; Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description; This function generates the linear code generators 
% the polynomial is of the form 1 + sum(x(i)) where 
% smn(x(i)) means that the output of the i'th cell of the register 
% is used as the input to the modulo-2 adder (exclusive or) 

% and the 1 means that the output of the adder is fed to the first cell 
% 

% so 1 + xl + x3 + xlO could be generated with poly =[1 13 10] and R 
% is the register to perform the operation on 

function r = gen_poly(poly, R) 

result = mod2( siun(R(poly(2;length(poly))))); % xor registers according to 

% polynomial 

r(2;length(R)) = R(1 ;(length(R)-l)); % shift registers right 

r(poly( 1 ))=result; % put result in register according 

% to polynomial 


% End of File gen_poly.m 


E. IN1T.M 


% File Name; initm 
% Last Updated; 11 Nov 96 
% Written By; T. H. Newman 
% Operating Environment: Sun Openwin & Win 95 
% Matlab Version 4.2c. 1 
% 

% Description: This script file sets the parameters for 
% the GPS signal simulation 

global w; 

% my_glb; 
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fo = 10.22999999543e6; 

fl = 120*fo; % frequency of L2 signal 

f2 = 154*fo; % frequency of LI signal 

To = .001; 

codejength =10; 

[s 1 .s2,s31=sig genfcode length); 
sigl=[(l:length(sl))',sr]; 
sig2=[(l :length(s2)y,s2']; 
sig3=[(1: Iength(s3))',s3']; 

kal_init; 

% End of File init.m 


F. KAL_INIT.M 


% File Name: kal_mit.m 
% Last Updated: 27 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Initialization file for Kalman Filter 
% 


global Pkkml; 
global R; 
global Q; 
global I; 
global Phi; 
global xkkml; 
global Tb; 
global Noise_gain; 
global dt; 

global w; 

% my_glb; 

A=[0 1; -w^2 0]; % xl = cos(wt) 

% x2 = xldot = -w*sin(wl) 
% x2_dot = -w^2*cos(wt) 
%x= [xl x2r 
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B=[0 w^2]’; 


% x_dot = A*[xl_dot, x2_dot] 
% =>A=[0 1;-w'^2 0]’ 

% x_dot = Ax + Bu 
% u is control input 


C=[l 0]; % will use to get xl estimate 

D=0; 


I=e}^e(2); % Identity Matrix 

[Phi,Del]=c2d(A,B,dt); % converts the continuous-time system: 

% x_dot = Ax -I- Bu 
% to the discrete-time state-space 
% system: 

% x[n+l] = Phi * x[n] + Del * u[n] 
% assuming a zero-order hold on tlie 
% inputs and sample time dt 


Pkkml=le6*I; 


R=l; 


%Q=0*I; 


% Initialize P(k|k-1) =M 

% where M = covariance of initial state 
% want this number large to give 
% weight to 1st measurment 
% var[e(0|0)] = var[v(0)] = R 
% V is reference input 
% (noise in measurement) 

% Q is measure of exitation by forcing function 


var_w=150; 

Q = [ (dtM)/4 (dt^3)/2, (dt'^3)/2 dt^2 ]*var_w; 


xkk = [0 0; 0 0]; % initialize x(k|k) to all zeros 

% 


xkkps= [0 0; 0 0]; 


x(k|k) is state of plant at time k given k observ ations 


xkkml = [0,0; 0,0]; 


% initialize x(k|k-l) => x(0|-l) = x_o = 0 


O/^- 

% clear all the variables used for saving plotting information 

clear myrand; 

clear x; 

clear sm; 

clear x_delay; 

clear z; 

clear y; 

clear xkkjjs; 

clear t; 
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% End of File kal init.m 


G. KALMAN.M 


% File Name: kalman.m 
% Last Updated: 28 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% QPSK signal 
% 

% Based on Professor Titus' algorithim which follows Prof Kirk's 
% NPS lecture notes " Optimal Estimation: An Introduction to the 
% Theory and Applications" 1975 All equation numbers are referenced 
% to the lecture notes. 

fimction y=kalman(z) 

global Pkkm 1; % Values initialized in kal_init 

global R; 

global Q; 

global I; 

global Phi; 

global xkkml; 

global w; 


C=ll,0]; 

Cps=[0,l]; 


G=Pkkml*C'*inv(C*Pkkml*C' +R); 


% Calculate Gain at time k 


% G(k) = P(k|k- 


l)*C'*inv(C*P(k|k-l)*C' + R) 


% eqn4.3-102a 


Pkk=(I-G*C)*Pklanl; % Calculate P(k|k) 


% P(k|k) = [I -G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 


Pkkml=Phi*Pkk*Phi' + Q; %P(k|k-l) 


eO = z-C*xkkml; 


xO=xkkml+G*eO; 


el = z-Cps*xkkml/w'; xl=Cps*(xkkml+G*el)/w; 


e2 = z+C*xkkml; 
e3 = z+Cps*xkkml/w; 


x2=-(xkkm 1+G*e2); 
x3= -Cps*(xkkml+G*e3)/w'; 


if abs(eO)<=abs(el); 
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e==eO; % Smallest error is with 

xkk = xkkml +G*(e); % old phase 

xkkps = xkk; 
else; 

e=el; % Smallest error is +90 

xkk = xkkml +G*(e); % is 
xklq>s = Cps*xkk/w; 
end 

if abs(e2) <= abs(e) 
e=e2; 

xkk = xkkml + G*(e); 

xkkps = -xkk; % Smallest error is -90 

end; 

if abs(e3) <= abs(e) 
e=e3; 

xkk = xkkml + G*e; % Smallest error is 180 

xkkps = -Cps*xkk/w; 
end; 


xkkml =Phi*xkk; % eqn 4.3-89a 

yl= C*xkk; % report the non shifted version 

)=[yl(l,l); xkkps(l); x0(l); xl(l); x2(l); x3(l)]; 


% End of File kalman.m 


H. KALMAN__U.M 


% File Name: kalman_u.m 
% Last Updated: 26 May 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% unbalanced QPSK signal 
% 

% Based on Professor Titus’ algorithim which follows Prof Kirk’s 
% NPS lecture notes ” Optimal Estimation: An Introduction to the 
% Theory and Applications” 1975 All equation numbers are referenced 
% to the lecture notes. 
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function y=kalman(z) 


global Pkkm 1; % Values initialized in kal_init 

global R; 

global Q; 

global I; 

global Phi; 

global xkknil; 

global w; 

global Al; 

global A2; 

C=[1,0]; 

Cps= [0,1]; 

G=Pkkml*C'*inv(C*Pkkml*C' +R); 
l)*C’*inv(C*P(k|k-l)*C' + R) 


Pkk=(I-G*C)*Pkkni 1; % Calculate P(k|k) 

% 

Pkkml=Phi*Pkk*Phi'+ Q; %P(k|k-l) 


% Calculate Gain at time k 

% G(k) = P(k|k- 

% eqn 4.3-102a 

P(k|k) = [1 -G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 


F=( 1/A2); % Factor for imbalanced quadature phase term 

cos_\it=C*xkkml; 

sin_wt=-Cps*xkkml/w; 


eO = z-(cos_wt + sqrt(2)*sin_wt); 
el = z-(cos_wt - sqrt(2)*sin_wt); 
e2 = z+(cos_wt - sqrt(2)*sin_wt); 
e3 = z+(cos_wt + sqrt(2)*sin_w1); 

if abs(eO)<=abs(el); 

e=e0; % Smallest error is with 

xkk = xkkml+G*(e); % old phase 

xkkps = (cos_wt+sqrt(2)*sin_wt); 
else; 

e=el; % Smallest error is +90 

xkk = xkkml +G*(e); % is 
xkkps = (coswt - sqrt(2)*sin_wt); 
end 


if abs(e2) <= abs(e) 
e=e2; 

xkk = xkkml + G*(e); 

xkkps = (cos wt - sqrt(2)*sin_wt); % Smallest error is -90 

end; 

if abs(e3) <= abs(e) 
e=e3; 
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xkk = xkkml + G*e; % Smallest error is 180 

xkkps = (cos_wt + sqrt(2)*sin_wt); 


end; 


xkkml=Phi*xkk; 


% eqn 4.3-89a 


yl= C*xkk; 


% report the non shifted version 


>-=[yl(l,l); xkkps(l)]; 


% End of File kalman_u.m 


L KALMAN2.M 


% File Name: kalman2.m 
% Last Updated: 3 May 97 
% Written B>’: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% BPSK signal. 

% 

% Based on Professor Titus' algorithim which follows Prof Kirk’s 
% NPS lecture notes" Optimal Estimation: An Introduction to the 
% Theory and Applications" 1975 All equation numbers are referenced 
% to the lecture notes. 

function y=kalman(z) 

global Pkkml; % Values initialized in kal_init 

global R; 

global Q; 

global 1; 

global Phi; 

global xkkml; 

global w; 


C=[1,0]; 
Cps =10,11; 


G=Pkkml*C'*inv(C*Pkkml*C’ +R); 


% Calculate Gain at time k 


% G(k) = P(k|k- 


l)*C’*inv(C*P(k|k-l)*C + R) 


% eqn4.3-102a 
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Pkk=(I-G*C)*Pkknil; % Calculate P(k|k) 

% P(k|k) = [I -G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 

Pkkml=Phi*Pkk*Phi' + Q; % P(k|k-1) 


eO = z-C*xkkml; xO=xkkml+G*(eO); 
e2 = z+C*xkkml; x2=-(xkknil+G*e2); 


if abs(eO)<= abs(e2); 

e=eO; % Smallest error is with 

xkk = xkkml +G*(e); % old phase 

xkkps = xkk; 
else 
e=e2; 

xkk = xkkml + G*(e); 

xkkps = -xkk; % Smallest error is -90 

end; 


xkkml=Phi*xkk; %eqn4.3-89a 

yl= C*xkk; % report the non shifted version 

y=( 5 'l(l,l); xkkps(l)]; 


% End of File kalman2.m 


J. MOD2.M 


% File Name; mod2.m 
% Last Updated: 18 Oct 96 
% Written B>': T. H. Newman 
% Operating Enviromnent: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Returns the modulo 2 value of input 
% if input is multiple of 2 retimis 0 
% else returns 1 
% 

% 


ftmction x = mod2(y) 


if (y/2)* 10 == round(y/2)* 10 
x = 0; 
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else 

x=l; 

end 


% End of File mod2.ni 


K. MPATH.M 


% File Name: mpath.m 
% Last Updated: 28 Apr 97 
% Written By: T. H. Newman 
% Operating En\iromnent: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function simulates the multipath environment 

function sm = mpath(s) 

global delay 
global alpha 
global tau 

len=length(delay); 

if (len > 1) 

delay (2:len) = delay(l :len-l); % Shift delay vector 
delay(l) = s; 

sm = sum(delay(tau).*alpha); 
else 

delay(l)=s; 

sm=s; 

end; 


% End of File mpath.m 


L. MPATH_AP.M 


% File Name: mpath_ap.m 
% Last Updated: 4 Jun 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 
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% Description; Kalman filter for multipath %vith a priori 
% information. 

clear; % clear old variables 

w=2; % = 2*pi*f 

Tf=15; % finish time ~ sec 
dt=.l; 

noise_fac=0; % change scale on noise 

%multipath a priori information 

alpha = .5; % delayed signal amplitude 

tau = 10 ; % time of delay 

Amp = 1; % amphtude of line of sight signal 


A=[0 1 0 0; -w^2 0 0 0;0 0 0 1 ;0 0 -w^2 0]; 
B=[0 ^"^2 0 0]'; 

C=[l 0 1 0]; 

[Phi,Del]=c2d(A,B,dt); 

I=eye(4); 

Pkkml=le6*I; 

R=l; 

Q=0.1*I; 
kmax=Tf/dt +1 ; 


switch 1 = 3; % time of phase changes 

svAitch2 = 9; 

time=0; 

for (i=l ;kmax) % generate line of sight signal 

t(i)=time; 
if t(i) <= switchl 

los(i) = cos(w*t(i)); % line of sight direct path signal 
elseif t(i) >= switchl & t(i) <= switch2 
los(i) = sin(w*t(i)); 
else 

los(i) = cos(w*t(i)); 
end; 

time=time+dt; 
end; 
tl=t; 
time = t; 
clear t; 

for i= 1 ;length(los)-tau % generate delayed signal 

x(i) = los(i+tau); % Line of sight signal x(t) 

delay(i) = los(i)*alpha; % delayed signal alpha*x(t-tau) 
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z(i) = x(i) + dela 3 '(i); % z = x(t) + alpha*x(t-tau) 
t(i) = time(i+tau);; 
end; 

z = z+noise_fac*randn(size(z)); % add in noise 

wt = w*time(tau+l); 
wtmtau = w*time(i); 
xkkml=[ cos(wt) 

-w*sin(w1) 

alpha*cos(wtmtau) 

-w*alpha*sin( wtmtau)]; 


for (i=l;length(z)) 

G=Pkkml*C*inv(C*Pkkml*C +R); 
Pkk=a-G*C)*Pkkml; 
Pkkml=Phi*Pkk*Phi' + Q; 
eO=z(l,i)-[l 0 1 0]*xkkml(;,i); 

el=z(l,i)-[0 -1/w 1 0]*xkkml(:,i); 
e2=z(l,i)-[0 -1/w 0 -l/w]*xkkml(;,i); 
e3=z(l,i)-[i 0 0 -l/w]*xkkml(:,i); 
%e_vect(:,i) = [ eO, el, e2, e3 ]'; 


if abs(eO)<=abs(el) & abs(eO) <= abs(e2) % abs(e0)<=abs(e3); % t(i)<=sl; 

%if (t(i) <= s\\itchl) 
e=eO; 

xkk(:,i) = xkkml(:,i) +G*(e); 
xkkps(l,i) = xkk(l,i); 
xkkps(3,i) = xkk(3,i); 

elseif abs(el) <= abs(eO) & abs(el) <= abs(e2) & abs(el)<=abs(e3) % si < t < sl+phi 
%elseif (switchl <= t(i) & t(i) <= switchl+tau*dt) 
e=el; 

xkk(: ,i)=xkkm 1 (:,i)+G*(e); 
xkkps(l,i) = -xkk(2,i)/w; 
xkkps(3,i) = xkk(3,i); 

elseif abs(e2)<=abs(e0) & abs(e2)<=abs(el) & abs(e2)<=abs(e3) % sl+tau < t < s2 
%elseif (switchl+tau*dt <= t(i) & t(i) <= sv\itch2) 
e=e2; 

xkk(:,i)=xkkml (:,i)+G*e; 
xkkps(l,i) = -xkk(2,i)/w; 
xkkps(3,i) = -xkk(4,i)/w; 

elseif abs(e3)<=abs(eO) & abs(e3)<=abs(el) & abs(e3)<=abs(e2) 

%elseif (switch2 <= t(i) & t(i) <= switch2+tau*dt) 

% s2 < t < s2+phi 

e=e3; 

xkk(; ,i)=xkkml (;,i)+G*e; 
xkkps( 1 ,i)=xkk( 1 ,i); 
xkkps(3 ,i)=-xkk(4,i)/w; 
else 

e=eO; 
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xkk(:,i) = xkkml(:,i) +G*(e); 
xldq}s(l,i) = xkk(l,i); 
xkkps(3,i) = xkk(3,i); 

end; 

xkkml(;,i+l)=Phi*xkk(;,i); 

end; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%% 


figure; 

mpath=delay; 

los=x; 

figure; 

subplot(3,l,l); 

plot(t,los,’r’); 

ylabel('los’); 

subplot(3,l,2); 

plot(t,mpath,'b'); 

ylabel('mpath'); 

subplot(3,l,3); 

plot(t,z,'g'); 

ylabel('z'); 

figure; 

subplot(2,l,l) 

plot(t,los,’r’,t,xkkps(l,;),'b-.’); 
title('los & xklq3s(l,:) (estimate = 

subplot(2,l,2); 

plot(t,mpath,'r',t,xkkps(3,;), !>.'); 
title('mpath & xkkps(3,;)'); 


% End of File mpath ap.m 
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M. 


MY GLB.M 


% File Name: my_glb.m 
% Last Upxiated: 16 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Global declaration for gps simulation and 
% kalman filter. 


global C; % Values initialized in kal_init.m 

global Pkkml; 
global R; 
global Q; 
global I; 
global Phi; 
global xkkml; 

global w; 

n = i; 

To=l; 
w=2*pi*fl; 
phase_init = pi/6; 
num_save = 100000; 


% NOTE: These are the same as in init.m 
% and are placed here for use in non-gps model. 

% Set frequency for simulation 
% Phase offset for sine wave in BPSK modulator 
% Number of points to save to workspace for plotting later 


% End of File my_glb.m 


N. SHIFT_REG.M 


% File Name: shift_reg.m 
% Last Updated: 21 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This fimction generates the shift register 
% linear code generators where R = the input values 
% size = register size. The register is shifted right 
% and returned in a new vector 

fimction r = shift_reg(R, size) 

r = ones(l,size) 


% End of File shift_reg.m 
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o. 


SIG GEN.M 


% File Name: si g g en.m 
% Last Updated: 11 Nov 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description; Returns either 1 P(Y) code xor with data 
% 2 P(Y) code 

% or 3 C/A code xor with data 

% the number is the sig_num which determines which signal 
% is returned 
% 


function [sl,s2,s3] = sigjgen(code_length) 
sv=l; 

xlepoch = 0; 
xlA_epoch = 0; 
xlB_epoch = 0; 
x2A_epoch = 0; 
x2B_epoch = 0; 

%code_length = 10; 

r_gl = ones( size( 1:10));% initialize register values 

r_g2 = ones( size( 1:10)); 

r_xla =[000 100100100]; 

r_xlb =[001010101010]; 

r_x2a=[10 100 1 00 100 1]; 

r_x2b =[001010101010]; 

shift_register = zeros(size(l:37)); 

t = 0; 

dataclock = 20; 


for i = 1 :code_length 
if data_clock == 20 

data = 1; % set to all ones to ignore 

data_clock = 0; 
else 


end; 


data_clock = data_clock + 1; 


forj =1:10 

tr_gl, g_ll = gl(r_gl, xl_epoch); 


% gen C/A code at 1.023 MHz 
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[r_g2, g_2] = g2(r_g2, xl_epoch, sv); 
g = xor(g_l, g_2); 

for k=l:10 % gen P/A code at 10.23 MHz 

t=t+l; % or 10 * C/A code generation rate 

[r_xla, A] = xla(r_xla, xlA_epoch); 

[r xlb, B] = xlb(r_xlb, xlB_epoch); 

[r_x2a, C] = x2a(r_x2a, x2A_epoch); 

[r_x2b, D] = x2b{r_x2b, x2B_epoch); 

xl=xor(A,B); 

x2=xor(C,D); 

shift_register(2;37) = shift_register(l:36); 
shifl_register(l) = x2; 
p=xor(xl, shift_register(sv)); 

sl(t) = xor(p,data); 
s2(t) = p; 

s3(t) = xor(g, data); 
end; % for k; 

end; %forj 

end; % for i 


% End of File sig gen.m 


P. TEST_CA.M 


% File Name: test_CA.m 
% Last Updated: 31 Mar 97 
% Written B>': T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to test C/A code generation. Output 
% should match First 10 C/A-Chips (Octal) of GPS-ICD-2000 table 

for sv=l:37 

r_gl = roimd( rand( size(l: 10))); % initialize to some random state 
r_g2 = round( rand( size(l:10))); 
xl = 1; 
i=l; 

[r_gl, g_l(i)] = gl(r_gl, xl); % Initialize at epoch 

[r_s2, g_2(i)] = g2(r_g2, xl, sv); 
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% Simulate 1st 10 chips 


xl=0; 
for i=2:10 

[r_gl. g_l(i)] = gl(r_gl, xl); 
[r_g2, g_2(i)] = g2(r_g2, xl, sv); 
end; 


g = xor(g_l, g_2); 


% 


digl = g(l); 

dig2 = g(2)*4 + g(3)*2 + g(4); % convert to octal number 

dig3 = g(5)*4 + g(6)*2 + g(7); 

dig4 = g(8)*4 + g(9)*2 + g(10); 

result(sv) = digl*10'^3 + dig2*10'^2 + dig3*10 + dig4; 

end 

result % display results 

% End of File test CA.m 


Q. TEST_P.M 


% File Name: test_P.m 
% Last Updated: 1 April 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to test P code generation. Output 
% should match First 12 P-Chips (Octal) of GPS-ICD-2000 table 
% Note that only the first 13 sv's are calculated. SV 9 - 37 are 
% identical 

for s\' =1:13 

r_xla = round( rand( size(l:12))); % initialize to some random values 
r_xlb= round( rand( size(l:12))); 
r_x2a = round( rand( size(l:12))); 
r_x2b = round( rand( size(l: 12))); 

epoch = 1; 
i=l; 

[r_xla, A(i)] = xla(r_xla, epoch); 

[r xlb, B(i)] = xlb(r_xlb, epoch); 

[r_x2a, C(i)] = x2a(r_x2a, epoch); 

[r_x2b, D(i)] = x2b(r_x2b, epoch); 


epoch=0; 

fori=2:12 

[r_xla, A(i)] = xla(r_xla, epoch); 
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[r xlb, B(i)] = xlb(r_xlb, epoch); 

[r_x2a, C(i)] = x2a(r_x2a, epoch); 

[r_x2b, D(i)] = x2b(r_x2b, epoch); 
end; 

xl=xor(A,B); 

x2=xor(C,D); 

delay=sv; 

xl=[xl, zeros(l,delay)]; 
x2=[ones(l,delay), x2 ]; 
p=xor(xl, x2); 

digl=p{l)*4 + p(2)*2 + p(3); 
dig2 = p(4)*4 + p(5)*2 + p(6); 
dig3 = p(7)*4 + p(8)*2 + p(9); 
dig4 = p(10)*4 + p(ll)*2 +p(12); 
result(sv) = digl*10^3 + dig2*10^2 + dig3*10 + dig4; 
end; % for sv 
result 

% End of File test P.m 


R. in_iNrr.M 


% File Name: ui_init.m 
% Last Updated: 4 May 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Called by ui_vars.m to initialize parameters 
% after the Init button is pushed. 

% 


delay=zeros(size(r.max(tau))); 
ifSNR>0 

Noise_gain = sqrt((Amp^2)/( 10^(SNR/10))); 
else 

Noise_gain= 1; 
end; 

if SNR >=100 
Noise_gain = 0; 
end; 

w=2*pi*f; 

disp(sprintf(' f=%2.2f Initial Phase = %2.2f, f, p_degree)) 
disp(sprintf(' Amp=%2.2f, Amp)) 
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disp(sprintf(' dt=%2,2f, stop = %2.2f, Tb=%2.2f, dt, stop, Tb)) 
disp(sprintf(' SNR=%2.2f dB, Noise_gain=%2.2f,SNR, Nois e g ain)) 
disp('tau='), disp(tau) 
disp('alpha='),disp(alpha) 


% End of File ui init.m 


S. UI_VARS.M 


% File Name: ui_vars.m 
% Last Updated: 28 Apr 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: An Editable Text windows to control 
% initialization values for kal_init 
% Calls: my_global.m to initialize values 
% kal_init.m to initialize kalman filter variables 

% init.m to generate GPS code sample 

global MainFig 

% below are the initial values that will be in the text window 

global C; % Values initialized in kal_init.m 

global Pkkml; 

global R; 

global Q; 

global I; 

global Phi; 

global xkkml; 

global Tb; 

global Noise_gain; 

global dt; 

global Amp; % Amplitude for BPSK 

global tau; 

global alpha; 

global delay; 

global w; 

global clock; 

my_glb; 

% This initilizes w, phase init, num_save, Tb, Noise gain, Al, A2, dt 
% Retreive initial values from the strings 
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start = 0; %start time 

stop = 10; %stop time 

f = w/(2*pi); %signal frequency 

p degree = phase_mit* 180/pi; %phase in degrees 

SNR= 100; 


f_string = sprintfC %2.2f,f);; 
p_degree_string = sprintf(' %2.2f ,p_degree); 

Amp_string = sprintfC %2.2f ,Amp); 
dt_string = sprintfC %2.2f ,dt);; 

start_string = sprintfC %2.2f,start); 

stop_string = sprintfC %2.2f ,stop); 

Tb_string = sprintfC %2.2f,Tb); 

SNR_string = sprintfC %2.2f, SNR); 


tau_string = ' [1 ]'; 

alpha_string =' [1]'; 

tau=eval(tau_string); 
alpha = eval(alpha_string); 

%define the Main Figure Window. 

num_entries = 9; 

sep=45; 

ymax = sep*(num_entries)+10; 

x_pos=25; y_pos=35; dx=125; dy = ymax; 

MainFig = figure ('Position', [x_pos y_pos dx 
'Color', 'white','NmnberTitle','off,'Name',... 
'Simulation Variables',... 

'MenuBar','none'); 


0 /„- 

% Define UIMENU for pull down menu options 

e=uimenuClaber,'lnit','callback',... 

'kal_init,ui_init;'); 

f=uimenuClaber, 'Plot'); 

uimenu(f, 'laber,'Phase Lock Loop','callback','plot_k'); 

uimenu(f,'laber,'Delay','callback','plot_d'); 

uimenu(f,'laber,'QPSK','callback','plot_q2'); 

uimenu(f,'label','BPSK','caUback','plotbpsk'); 

uimenu(f,'laber,'Error','callback','plot_err'); 

uimenu(f,'laber,'Variables','callback','plot_k2'); 

g=uimenuClaber, 'Sim'); 
gl=uimenu(g, 'label','BPSK'); 
uimenu(gl, 'label', 'BPSK Simulation','callback','kal_sim'); 
uimenu(gl, 'label', 'PLL Simulation','callback','pH'); 
g2=uimenu(g, 'label', 'QPSK'); 
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uimenu(g2, 'label','QPSK Simulation’, 'callback','kalq_sim’); 
uimenu(g2, 'label', 'PLL Simulation','callback’,'pll_q'); 
g3=uimenu(g, 'label', 'Unbalanced QPSK'); 
uimenu(g3, 'laber,'Unbalanced QPSK Simulation', 'callback',’kalu_sim'); 
uimenu(g3, 'label', 'PLL SimuIation','callback','pll_u'); 
g4=uimenu(g,'label','Multipath','callback','mpathsim'); 


n=l; %uimenu number 

jinax = dy; 

xpos=10; 

dx = 100; 

dy=20; 

space = 22; 

label_fg_color = 'Black'; 
label_bg_color = 'White'; 
box fg color = 'Black'; 
box_bg_color = 'Yellow'; 

f_label = uicontrol(MainFig, ’St)'le',’text', 'Pos',[xpos ymax-n*seFH-space dx cfy],... 
'String', 'Signal freq (Hz)','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 


f_box = uicontrol(MainFig,'Style','edit','String', f_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

Toregroundcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

’CaUback',’f = eval( get(f_box,"String")); w-2*pi*f;'); 

n=n+l; 

phase_label = uicontrol(MainFig, 'Style',text', 'Pos',[xpos ymax-n*seiri-space dx dy],... 
'String', 'Phase Offset °','ForegroundColor',label_fg_color,... 
'BackgroimdColor',label_bg_color); 

phase_box = uicontrol(MainFig,’Style','edit','String', p_degree_string,... 

'Pos’,[xpos ymax-n*sep dx cfy],... 

'Foregroxmdcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

'Callback','p_degree = eval( get(phase_box,"String") );phase_init=p_degree*pi/180'); 
n=n+l; 

amp label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos ymax-n*seiri-space dx dy],... 
'String', 'Amplitude','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

amp box = uicontrol(MainFig,’Style','edit','String', Amp_string,... 

'Pos',[xpos ymax-n*sep dx d\'],... 
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'Foregroundcolor',box_fg_color,... 
'BackgroundColor',box_bg_color,... 

'Callback','A1 = eval( get(amp_box,"String"));'); 


n=n+l; 

samp label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos yinax-n*sep+space dx cfy], 
'String', 'Ts (sec)',ToregroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

sampbox = uicontrol(MainFig,'Style','edit','String', dt_string,... 

Pos',[xpos ymax-n*sep dx dy],... 

T oregroundcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

'Callback','dt = eval( get(samp_box,"String") 

n=n+l; 

stopjabel = uicontrol(MainFig, 'Style','text', Pos',[xposyinax-n*sep+space dx d)'],.. 
'String', 'Stop time (sec)','ForegroimdColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

stop_box = uicontrol(MainFig,'Style','edit','String', stop_string,... 

'Pos',[xpos yniax-n*sep dx dj'],... 

'Foregroundcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

'Callback','stop = eval (get(stop_box,"String"));'); 

n=n+l; 

Tb_label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos ymax-n*sep+space dx dy],... 
'String', 'Tb (sec)','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

Tb_box = uicontrol(MainFig,'Style','edit','String', Tb_string,... 

'Pos',[xpos ymax-n*sep dx (fy],... 

'Foregroundcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

•Callback',’Tb = eval (get(Tb_box,"String"));'); 


n=n+l; 

SNR_label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos }Tnax-n*sep+space dx dy],. 
'String', 'SNR dB','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

SNR_box = mcontrol(MainFig,'Style','edit','String', SNR_string,... 

'Pos',[xpos ymax-n*sep dx c^],... 

'Foregroundcolor',box_fg_color,... 

'BackgroimdColor',box_bg_color,... 

'Callback','SNR = eval( get(SNR_box,"String"));'); 
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n=n+l; 


tau_label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos yinax-n*sep+space dx dy],... 
'String', 'tau','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg^color); 

tau_box = uicontrol(MainFig, 'Style','edit','String', tau string,... 

'Pos',[xpos ymax-n*sep dx <ty],... 

'Foregroimdcolor',box_fg_color,... 

'BackgroimdColor',box_bg_color,... 

'Callback','tau = eval( get(tau_box,"String"));'); 

n=n+l; 


alphajabel = uicontrol(MainFig, 'Style','text', 'Pos',[xpos ymax-n*sep+space dx dy],... 
'String', 'alpha','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

alphajbox = uicontrol(MainFig,'Style','edit','String', alpha_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

'Foregroundcolor',box_fg_color,... 

'BackgroundColor',box_bg_color,... 

'Callback','alpha = eval( get(alpha_box,"String"));'); 

n=n+l; 

0 /„ - 0 ^ 

% Define the UICONTROL button for callbacks. 

% bl=uicontrol('style','pushbutton','string','Initialize',... 

% 'position', [30 ymax-n*sep 65 25],'Callback',... 

% 'kal_init; delay=zeros(size(l:max(tau)));'); 


% End of File ui vars.m 


T. X1A.M 


% File Name: xla.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 
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% Description: This function generates the xla(t) signal 
% for GPS P Code generator 

function [new R, A] = xla(R, xla_epoch) 

if xla_epoch==l 

R = [0 0010010010 0]; % Reset at XIA epoch 
end; 

new_R = genjx)ly([l,6,8,ll,12], R); 

A=R(12); 

% End of File xla.m 


U. X1B.M 


% File Name: xlb.m 
% Last Updated: 18 Oct 96 
% Written By; T. H. Newman 
% Operating Environment: Window's 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the xlb(t) signal 
% for GPS P Code generator 

function [new'_R, B] = xlb(R, xlb_epoch) 

if xlb_epoch=l 

R= [00 1 0 1 0 1 0 1 0 1 0]; % Reset at XIB epoch 
end; 

new_R = genjx)ly([l 1 2 5 8 9 10 11 12], R); 

B=R(12); 

% End of File xlb.m 


V. X2A.M 


% File Name: x2a.m 
% Last Updated; 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the x2a(t) signal 
% for GPS P/A Code generator 
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function (new_R, C] = x2a(R, x2a_epoch) 


if x2a_epoch==l 

R = [10100100100 1]; % Reset at X2a epoch 
end; 

new_R = gen_poly([l 1 3 4 5 7 8 9 10 11 12], R); 
C=R(12); 

% End of File x2a.m 


W. X2B.M 


% File Name; x2b.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment; Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the x2b(t) signal 
% for GPS P Code generator 

function [new_R, D] = x2b(R, x2b_epoch) 


if x2b_epoch==l 

R= [00 1 0 1 0 1 0 1 0 1 0]; % Reset at X2B epoch 
end; 

new_R = genj30ly([l 2 3 4 8 9 12], R); 

D=R(12); 

% End of File x2b.m 
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APPENDIX B SEMULINK DIAGRAMS 
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E. KALQ_SIM.M 
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APPENDIX C EXAMPLE SNR CALCULATION 


X 

dB2R(x) = 10^*^ convert dB to ratio dB(x) =101og(x) convert ratio to dB 


Minimum Received Power 

L1 with C/A Code 

C:=-160.0 ~dB_W 

Boltzman's Constant 

k =-228.6 -dBW/Hz-K 

Noise Temperature 

Tn :=28 ~ dB K 

Received Noise 

No :=k-i- Tn 

No =-200.6 

Carrier to Noise Density Ratio 

CNR:=C- No 

CNR=40.6 ~dbHz 


Chip Data Rate 

R = 1.02310^ ~ chips/sec 

Information Data Rate 

R_db :=dB(R) 

R_db =60.099 

R_i :=50 ~ bits/sec 

Band Width - QPSK 

B;=0.6R 

B_dB =dB(B) 

B_dB = 57.88 -dBbps 

Signal to Noise Ratio 

SNR = CNR- B_dB 

SNR = -17.28 -dB 

Processing Gain 

PG:=dBi-^l PG=43.109 - dB 

\R_i/ 

Signal to Noise Ratio 

SNR :=SNR-hPG SNR = 25.829 - dB 
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